home *** CD-ROM | disk | FTP | other *** search
/ Monster Media 1996 #15 / Monster Media Number 15 (Monster Media)(July 1996).ISO / prog_c / cuj0696.zip / DWYER.ZIP / LIB / NDTRI.C < prev    next >
C/C++ Source or Header  |  1989-03-21  |  10KB  |  411 lines
  1. /*                            ndtri.c
  2.  *
  3.  *    Inverse of Normal distribution function
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * double x, y, ndtri();
  10.  *
  11.  * x = ndtri( y );
  12.  *
  13.  *
  14.  *
  15.  * DESCRIPTION:
  16.  *
  17.  * Returns the argument, x, for which the area under the
  18.  * Gaussian probability density function (integrated from
  19.  * minus infinity to x) is equal to y.
  20.  *
  21.  *
  22.  * For small arguments 0 < y < exp(-2), the program computes
  23.  * z = sqrt( -2.0 * log(y) );  then the approximation is
  24.  * x = z - log(z)/z  - (1/z) P(1/z) / Q(1/z).
  25.  * There are two rational functions P/Q, one for 0 < y < exp(-32)
  26.  * and the other for y up to exp(-2).  For larger arguments,
  27.  * w = y - 0.5, and  x/sqrt(2pi) = w + w**3 R(w**2)/S(w**2)).
  28.  *
  29.  *
  30.  * ACCURACY:
  31.  *
  32.  *                      Relative error:
  33.  * arithmetic   domain        # trials      peak         rms
  34.  *    DEC      0.125, 1         5500       9.5e-17     2.1e-17
  35.  *    DEC      6e-39, 0.135     3500       5.7e-17     1.3e-17
  36.  *    IEEE     0.125, 1        20000       7.2e-16     1.3e-16
  37.  *    IEEE     3e-308, 0.135   50000       4.6e-16     9.8e-17
  38.  *
  39.  *
  40.  * ERROR MESSAGES:
  41.  *
  42.  *   message         condition    value returned
  43.  * ndtri domain       x <= 0        -MAXNUM
  44.  * ndtri domain       x >= 1         MAXNUM
  45.  *
  46.  */
  47.  
  48.  
  49. /*
  50. Cephes Math Library Release 2.1:  January, 1989
  51. Copyright 1984, 1987, 1989 by Stephen L. Moshier
  52. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  53. */
  54.  
  55. #include "mconf.h"
  56. extern double MAXNUM;
  57.  
  58. #ifdef UNK
  59. /* sqrt(2pi) */
  60. static double s2pi = 2.50662827463100050242E0;
  61. #endif
  62.  
  63. #ifdef DEC
  64. static short s2p[] = {0040440,0066230,0177661,0034055};
  65. #define s2pi *(double *)s2p
  66. #endif
  67.  
  68. #ifdef IBMPC
  69. static short s2p[] = {0x2706,0x1ff6,0x0d93,0x4004};
  70. #define s2pi *(double *)s2p
  71. #endif
  72.  
  73. #ifdef MIEEE
  74. static short s2p[] = {
  75. 0x4004,0x0d93,0x1ff6,0x2706
  76. };
  77. #define s2pi *(double *)s2p
  78. #endif
  79.  
  80. /* approximation for 0 <= |y - 0.5| <= 3/8 */
  81. #ifdef UNK
  82. static double P0[5] = {
  83. -5.99633501014107895267E1,
  84.  9.80010754185999661536E1,
  85. -5.66762857469070293439E1,
  86.  1.39312609387279679503E1,
  87. -1.23916583867381258016E0,
  88. };
  89. static double Q0[8] = {
  90. /* 1.00000000000000000000E0,*/
  91.  1.95448858338141759834E0,
  92.  4.67627912898881538453E0,
  93.  8.63602421390890590575E1,
  94. -2.25462687854119370527E2,
  95.  2.00260212380060660359E2,
  96. -8.20372256168333339912E1,
  97.  1.59056225126211695515E1,
  98. -1.18331621121330003142E0,
  99. };
  100. #endif
  101. #ifdef DEC
  102. static short P0[20] = {
  103. 0141557,0155170,0071360,0120550,
  104. 0041704,0000214,0172417,0067307,
  105. 0141542,0132204,0040066,0156723,
  106. 0041136,0163161,0157276,0007747,
  107. 0140236,0116374,0073666,0051764,
  108. };
  109. static short Q0[32] = {
  110. /*0040200,0000000,0000000,0000000,*/
  111. 0040372,0026256,0110403,0123707,
  112. 0040625,0122024,0020277,0026661,
  113. 0041654,0134161,0124134,0007244,
  114. 0142141,0073162,0133021,0131371,
  115. 0042110,0041235,0043516,0057767,
  116. 0141644,0011417,0036155,0137305,
  117. 0041176,0076556,0004043,0125430,
  118. 0140227,0073347,0152776,0067251,
  119. };
  120. #endif
  121. #ifdef IBMPC
  122. static short P0[20] = {
  123. 0x142d,0x0e5e,0xfb4f,0xc04d,
  124. 0xedd9,0x9ea1,0x8011,0x4058,
  125. 0xdbba,0x8806,0x5690,0xc04c,
  126. 0xc1fd,0x3bd7,0xdcce,0x402b,
  127. 0xca7e,0x8ef6,0xd39f,0xbff3,
  128. };
  129. static short Q0[36] = {
  130. /*0x0000,0x0000,0x0000,0x3ff0,*/
  131. 0x74f9,0xd220,0x4595,0x3fff,
  132. 0xe5b6,0x8417,0xb482,0x4012,
  133. 0x81d4,0x350b,0x970e,0x4055,
  134. 0x365f,0x56c2,0x2ece,0xc06c,
  135. 0xcbff,0xa8e9,0x0853,0x4069,
  136. 0xb7d9,0xe78d,0x8261,0xc054,
  137. 0x7563,0xc104,0xcfad,0x402f,
  138. 0xcdd5,0xfabf,0xeedc,0xbff2,
  139. };
  140. #endif
  141. #ifdef MIEEE
  142. static short P0[20] = {
  143. 0xc04d,0xfb4f,0x0e5e,0x142d,
  144. 0x4058,0x8011,0x9ea1,0xedd9,
  145. 0xc04c,0x5690,0x8806,0xdbba,
  146. 0x402b,0xdcce,0x3bd7,0xc1fd,
  147. 0xbff3,0xd39f,0x8ef6,0xca7e,
  148. };
  149. static short Q0[32] = {
  150. /*0x3ff0,0x0000,0x0000,0x0000,*/
  151. 0x3fff,0x4595,0xd220,0x74f9,
  152. 0x4012,0xb482,0x8417,0xe5b6,
  153. 0x4055,0x970e,0x350b,0x81d4,
  154. 0xc06c,0x2ece,0x56c2,0x365f,
  155. 0x4069,0x0853,0xa8e9,0xcbff,
  156. 0xc054,0x8261,0xe78d,0xb7d9,
  157. 0x402f,0xcfad,0xc104,0x7563,
  158. 0xbff2,0xeedc,0xfabf,0xcdd5,
  159. };
  160. #endif
  161.  
  162.  
  163. /* Approximation for interval z = sqrt(-2 log y ) between 2 and 8
  164.  * i.e., y between exp(-2) = .135 and exp(-32) = 1.27e-14.
  165.  */
  166. #ifdef UNK
  167. static double P1[9] = {
  168.  4.05544892305962419923E0,
  169.  3.15251094599893866154E1,
  170.  5.71628192246421288162E1,
  171.  4.40805073893200834700E1,
  172.  1.46849561928858024014E1,
  173.  2.18663306850790267539E0,
  174. -1.40256079171354495875E-1,
  175. -3.50424626827848203418E-2,
  176. -8.57456785154685413611E-4,
  177. };
  178. static double Q1[8] = {
  179. /*  1.00000000000000000000E0,*/
  180.  1.57799883256466749731E1,
  181.  4.53907635128879210584E1,
  182.  4.13172038254672030440E1,
  183.  1.50425385692907503408E1,
  184.  2.50464946208309415979E0,
  185. -1.42182922854787788574E-1,
  186. -3.80806407691578277194E-2,
  187. -9.33259480895457427372E-4,
  188. };
  189. #endif
  190. #ifdef DEC
  191. static short P1[36] = {
  192. 0040601,0143074,0150744,0073326,
  193. 0041374,0031554,0113253,0146016,
  194. 0041544,0123272,0012463,0176771,
  195. 0041460,0051160,0103560,0156511,
  196. 0041152,0172624,0117772,0030755,
  197. 0040413,0170713,0151545,0176413,
  198. 0137417,0117512,0022154,0131671,
  199. 0137017,0104257,0071432,0007072,
  200. 0135540,0143363,0063137,0036166,
  201. };
  202. static short Q1[32] = {
  203. /*0040200,0000000,0000000,0000000,*/
  204. 0041174,0075325,0004736,0120326,
  205. 0041465,0110044,0047561,0045567,
  206. 0041445,0042321,0012142,0030340,
  207. 0041160,0127074,0166076,0141051,
  208. 0040440,0046055,0040745,0150400,
  209. 0137421,0114146,0067330,0010621,
  210. 0137033,0175162,0025555,0114351,
  211. 0135564,0122773,0145750,0030357,
  212. };
  213. #endif
  214. #ifdef IBMPC
  215. static short P1[36] = {
  216. 0x8edb,0x9a3c,0x38c7,0x4010,
  217. 0x7982,0x92d5,0x866d,0x403f,
  218. 0x7fbf,0x42a6,0x94d7,0x404c,
  219. 0x1ba9,0x10ee,0x0a4e,0x4046,
  220. 0x463e,0x93ff,0x5eb2,0x402d,
  221. 0xbfa1,0x7a6c,0x7e39,0x4001,
  222. 0x9677,0x448d,0xf3e9,0xbfc1,
  223. 0x41c7,0xee63,0xf115,0xbfa1,
  224. 0xe78f,0x6ccb,0x18de,0xbf4c,
  225. };
  226. static short Q1[32] = {
  227. /*0x0000,0x0000,0x0000,0x3ff0,*/
  228. 0xd41b,0xa13b,0x8f5a,0x402f,
  229. 0x296f,0x89ee,0xb204,0x4046,
  230. 0x461c,0x228c,0xa89a,0x4044,
  231. 0xd845,0x9d87,0x15c7,0x402e,
  232. 0xba20,0xa83c,0x0985,0x4004,
  233. 0x0232,0xcddb,0x330c,0xbfc2,
  234. 0xb31d,0x456d,0x7f4e,0xbfa3,
  235. 0x061e,0x797d,0x94bf,0xbf4e,
  236. };
  237. #endif
  238. #ifdef MIEEE
  239. static short P1[36] = {
  240. 0x4010,0x38c7,0x9a3c,0x8edb,
  241. 0x403f,0x866d,0x92d5,0x7982,
  242. 0x404c,0x94d7,0x42a6,0x7fbf,
  243. 0x4046,0x0a4e,0x10ee,0x1ba9,
  244. 0x402d,0x5eb2,0x93ff,0x463e,
  245. 0x4001,0x7e39,0x7a6c,0xbfa1,
  246. 0xbfc1,0xf3e9,0x448d,0x9677,
  247. 0xbfa1,0xf115,0xee63,0x41c7,
  248. 0xbf4c,0x18de,0x6ccb,0xe78f,
  249. };
  250. static short Q1[32] = {
  251. /*0x3ff0,0x0000,0x0000,0x0000,*/
  252. 0x402f,0x8f5a,0xa13b,0xd41b,
  253. 0x4046,0xb204,0x89ee,0x296f,
  254. 0x4044,0xa89a,0x228c,0x461c,
  255. 0x402e,0x15c7,0x9d87,0xd845,
  256. 0x4004,0x0985,0xa83c,0xba20,
  257. 0xbfc2,0x330c,0xcddb,0x0232,
  258. 0xbfa3,0x7f4e,0x456d,0xb31d,
  259. 0xbf4e,0x94bf,0x797d,0x061e,
  260. };
  261. #endif
  262.  
  263. /* Approximation for interval z = sqrt(-2 log y ) between 8 and 64
  264.  * i.e., y between exp(-32) = 1.27e-14 and exp(-2048) = 3.67e-890.
  265.  */
  266.  
  267. #ifdef UNK
  268. static double P2[9] = {
  269.   3.23774891776946035970E0,
  270.   6.91522889068984211695E0,
  271.   3.93881025292474443415E0,
  272.   1.33303460815807542389E0,
  273.   2.01485389549179081538E-1,
  274.   1.23716634817820021358E-2,
  275.   3.01581553508235416007E-4,
  276.   2.65806974686737550832E-6,
  277.   6.23974539184983293730E-9,
  278. };
  279. static double Q2[8] = {
  280. /*  1.00000000000000000000E0,*/
  281.   6.02427039364742014255E0,
  282.   3.67983563856160859403E0,
  283.   1.37702099489081330271E0,
  284.   2.16236993594496635890E-1,
  285.   1.34204006088543189037E-2,
  286.   3.28014464682127739104E-4,
  287.   2.89247864745380683936E-6,
  288.   6.79019408009981274425E-9,
  289. };
  290. #endif
  291. #ifdef DEC
  292. static short P2[36] = {
  293. 0040517,0033507,0036236,0125641,
  294. 0040735,0044616,0014473,0140133,
  295. 0040574,0012567,0114535,0102541,
  296. 0040252,0120340,0143474,0150135,
  297. 0037516,0051057,0115361,0031211,
  298. 0036512,0131204,0101511,0125144,
  299. 0035236,0016627,0043160,0140216,
  300. 0033462,0060512,0060141,0010641,
  301. 0031326,0062541,0101304,0077706,
  302. };
  303. static short Q2[32] = {
  304. /*0040200,0000000,0000000,0000000,*/
  305. 0040700,0143322,0132137,0040501,
  306. 0040553,0101155,0053221,0140257,
  307. 0040260,0041071,0052573,0010004,
  308. 0037535,0066472,0177261,0162330,
  309. 0036533,0160475,0066666,0036132,
  310. 0035253,0174533,0027771,0044027,
  311. 0033502,0016147,0117666,0063671,
  312. 0031351,0047455,0141663,0054751,
  313. };
  314. #endif
  315. #ifdef IBMPC
  316. static short P2[36] = {
  317. 0xd574,0xe793,0xe6e8,0x4009,
  318. 0x780b,0xc327,0xa931,0x401b,
  319. 0xb0ac,0xf32b,0x82ae,0x400f,
  320. 0x9a0c,0x18e7,0x541c,0x3ff5,
  321. 0x2651,0xf35e,0xca45,0x3fc9,
  322. 0x354d,0x9069,0x5650,0x3f89,
  323. 0x1812,0xe8ce,0xc3b2,0x3f33,
  324. 0x2234,0x4c0c,0x4c29,0x3ec6,
  325. 0x8ff9,0x3058,0xccac,0x3e3a,
  326. };
  327. static short Q2[32] = {
  328. /*0x0000,0x0000,0x0000,0x3ff0,*/
  329. 0xe828,0x568b,0x18da,0x4018,
  330. 0x3816,0xaad2,0x704d,0x400d,
  331. 0x6200,0x2aaf,0x0847,0x3ff6,
  332. 0x3c9b,0x5fd6,0xada7,0x3fcb,
  333. 0xc78b,0xadb6,0x7c27,0x3f8b,
  334. 0x2903,0x65ff,0x7f2b,0x3f35,
  335. 0xccf7,0xf3f6,0x438c,0x3ec8,
  336. 0x6b3d,0xb876,0x29e5,0x3e3d,
  337. };
  338. #endif
  339. #ifdef MIEEE
  340. static short P2[36] = {
  341. 0x4009,0xe6e8,0xe793,0xd574,
  342. 0x401b,0xa931,0xc327,0x780b,
  343. 0x400f,0x82ae,0xf32b,0xb0ac,
  344. 0x3ff5,0x541c,0x18e7,0x9a0c,
  345. 0x3fc9,0xca45,0xf35e,0x2651,
  346. 0x3f89,0x5650,0x9069,0x354d,
  347. 0x3f33,0xc3b2,0xe8ce,0x1812,
  348. 0x3ec6,0x4c29,0x4c0c,0x2234,
  349. 0x3e3a,0xccac,0x3058,0x8ff9,
  350. };
  351. static short Q2[32] = {
  352. /*0x3ff0,0x0000,0x0000,0x0000,*/
  353. 0x4018,0x18da,0x568b,0xe828,
  354. 0x400d,0x704d,0xaad2,0x3816,
  355. 0x3ff6,0x0847,0x2aaf,0x6200,
  356. 0x3fcb,0xada7,0x5fd6,0x3c9b,
  357. 0x3f8b,0x7c27,0xadb6,0xc78b,
  358. 0x3f35,0x7f2b,0x65ff,0x2903,
  359. 0x3ec8,0x438c,0xf3f6,0xccf7,
  360. 0x3e3d,0x29e5,0xb876,0x6b3d,
  361. };
  362. #endif
  363.  
  364. double ndtri(y0)
  365. double y0;
  366. {
  367. double x, y, z, y2, x0, x1;
  368. int code;
  369. double polevl(), p1evl(), log(), sqrt();
  370.  
  371. if( y0 <= 0.0 )
  372.     {
  373.     mtherr( "ndtri", DOMAIN );
  374.     return( -MAXNUM );
  375.     }
  376. if( y0 >= 1.0 )
  377.     {
  378.     mtherr( "ndtri", DOMAIN );
  379.     return( MAXNUM );
  380.     }
  381. code = 1;
  382. y = y0;
  383. if( y > (1.0 - 0.13533528323661269189) ) /* 0.135... = exp(-2) */
  384.     {
  385.     y = 1.0 - y;
  386.     code = 0;
  387.     }
  388.  
  389. if( y > 0.13533528323661269189 )
  390.     {
  391.     y = y - 0.5;
  392.     y2 = y * y;
  393.     x = y + y * (y2 * polevl( y2, P0, 4)/p1evl( y2, Q0, 8 ));
  394.     x = x * s2pi; 
  395.     return(x);
  396.     }
  397.  
  398. x = sqrt( -2.0 * log(y) );
  399. x0 = x - log(x)/x;
  400.  
  401. z = 1.0/x;
  402. if( x < 8.0 ) /* y > exp(-32) = 1.2664165549e-14 */
  403.     x1 = z * polevl( z, P1, 8 )/p1evl( z, Q1, 8 );
  404. else
  405.     x1 = z * polevl( z, P2, 8 )/p1evl( z, Q2, 8 );
  406. x = x0 - x1;
  407. if( code != 0 )
  408.     x = -x;
  409. return( x );
  410. }
  411.